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SECTION 1 


INTRODUCTION 


The present report includes work performed during the second year of 
a NASA-Lewis grant to study the energy absorption mechanisms during crack 
propagation in metal matrix composites. The first-year work was reported 
in Reference [1]. This previous report contains a literature review of 
micromechonics analyses of unidirectional composites, and a discussion of 
the relation of these prior studies to the present problem. 

During the first year, an existing elastoplast ie , finite element 
analysis and associated computer program [2,3] was used to predict the 
response of a unidirectional boron/aluminum composite to axial loading. 

For tills purpose, a longitudinal section model was constructed. This model 
permitted the study of the influence, of a broken fiber on the load distri- 
bution in adjacent unbroken fibers one and two fiber spacings away . It 
also permitted the determination of the rate of reloading of the broken 
fiber away from the site of the break. The influence of plastic deforma- 
tion of the aluminum matrix on the stress distributions was of special 
interest. The addition of a crack initiation and propagation capability 
was initiated, but minor programming difficulties prevented results being 
presented in the first-year report. 

The goal of the second-year study was to complete the crack propaga- 
tion addition and generate detailed numerical results. Also, it was desired 
to construct larger longitudinal section models, to determine stress 
redistributions and influences of a broken fiber beyond the second adjacent 
fiber, and to determine the extent of influence, of boundary loading 
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conditions. In addition, the existing generalized plane strain analysis 
was to be converted to an axisymmetrlc analysis, to permit the study of 
the response of a single fiber embedded in a cylindrical sheath of matrix 
material. 

Those second-year goals were achieved, as summarized in the next section, 
and described in detail in the remainder of this report. 



SECTION 2 


SUMMARY 


The crack initiation and propagation capability became operational 
early in the second-year study. The entire analysis was then converted from 
a Sigma 7 computer system to the University's new CDC Cyber 730/760 computer 
which had just been Installed. The much greater capacity and speed of 
this new system made it practical to analyze much larger finite element 
grid arrays. Thus, new models were constructed, involving a greater number 
of fibers adjacent to the broken fiber, and greater lengths along the 
fiber axes. The results obtained using these larger models were then com- 
pared with previous results. The results of the first-year study were also 
extended beyond first failure, to analyze crack propagation behavior. 

In addition to the longitudinal section models, transverse section 
models were also constructed and analyzed. These included arrays of un- 
broken fibers, and also a single broken fiber surrounded by unbroken libers. 
Although axi ;1 loading was of primary interest, transverse loading was also 
studied. This permitted the comparison of results with other transverse 
loading results available in the literature [4-6], for verification purposes 
The crack propagation capability of the micromechanics analysis was 
found to perform very well, and is now considered to be fully operational. 

The conversion of the generalized plane strain mioromechanics analysis 
to an axlsymmetric. formulation proved to be more difficult than anticipated. 
The difficulty was not in the basic reformulation, but In the detailed 
modifications required in the associated computer program. Ultimately a 
second program was developed, as a more practical approach than attempting 
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to make the ax i symmetric formulation an option in the original program. 

The axisymmetrie analysis and associated computer program is now operational. 
Several example problems are presented in this report to demonstrate its 
capabilities . 

While minor improvements will undoubtedly be incorporated into each of 
the now operational computer programs during the third-yenrd effort, they 
are essentially complete. Attention will thus be focused on correlating 
predictions with available experimental data, and making parametric studies 
of the influences of various experimental variables, such as fiber volume, 
matrix properties, locations of fiber breaks, etc. Also, use will be made 
of a new three-dimensional finite element analysis recently completed as 
part of another study [7], to examine further the three-dimensional nature 
of stress states around broken fibers. This analysis will be useful directly 
in analyzing practical problems, and also for verification purposes in 
establishing the limits of applicability of the two-dimensional generalized 
plane strain and axisymmetrie analyses, which are more economical to utilize. 

It is anticipated that NASA-Lewis will be generating specialized and 
carefully controlled supporting experimental data also, primarily using 
single fiber specimens, to study fiber fracture and matrix deformations. 

The various analyses will be correlated with fhese experimental data. 


SECTION 3 


GENERALIZED PLANK STRAIN ANALYSIS METHOD 

The analysis formulation was presented in detail In the first-year 
report [1]. This has not changed, the governing equations and the flow 
chart which defines the operational features of the computer program modi- 
fied to implement the analysis, which were presented in References [1]» 
remain valid. Thus, only a brief summary need be given here. 

The primary analytical tool used in the present study has been the 
micromechanics finite element analysis program developed by Miller and 
Adams [2,3]. This analysis was developed to investigate the microstress 
state in unidirectional composite materials subjected to axial and trans- 
verse mechanical loads, thermal gradients, and dilatation.il stresses due to 
moisture absorption by polymeric matrix materials. Among the special 
features of thin prior University of Wyoming analysis are its ability to 
model the elastoplastic stress-strain response of the isotropic matrix 
material, and in concert with the determination of thermal and moisture 
dllatational stresses, the functional dependence of the matrix material 
properties on temperature and moisture content. In other words, the 
elastic or plastic properties of any matrix material finite element are 
automatically computed to reflect the state of stress and the environmental 
conditions of temperature and humidity. The adjustment of material prop- 
erties is incorporated with the incremental loading technique that is em- 
ployed in this program. Once the initial temperature, moisture content, and/ 
or elastic stress level for the continuum have been specified, additional 
loads, be they mechanical or environmental, are introduced in increments 
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small enough to permit close approximation of the nonlinear matrix material 
properties by small linear segments. A detailed description of this technique 
is presented in References [1,2,3]. 

The bulk of the elastoplastic formulation in the present analysis pro- 
gram stems from previous work done by Adams [4-6, 8, 9]. The development of 
a generalized plane strain formulation, incorporation of the Branca solution 
technique [10], hygrothertnal loading, and material properties dependence 
on temperature and moisture was the subject of Miller's Ph.D. research [2], 
while the addition of crack initiation and crack propagation capability 
follows the approach developed by Adams [4-6], The analysis incorporates 
standard finite element techniques (see, among many other similar sources, 
References [11,12]). In fact, the primary organization and flow of the 
original computer program closely followed the suggestions of Appendix A 
of Reference [12], This flow and organization has subsequently been rather 
severely altered to include crack propagation capability. 

The finite element used in this study is a modified version of the 
familiar constant strain or simplex triangle. For this element, a linear 
displacement field within each element is assumed, to arrive at a functional 
representation of the potential energy of the system. The constant strain 
triangular element has some well-known limitations, but for the purposes of 
micromechanics analyses, it has proven to be an accurate, economical, and 
versatile tool. The trade-offs involved in the choice of the constant strain 
triangular element instead of one of the higher order finite elements is 
covered quite well by Miller and Adams in Chapter 3 of Reference [2], 

A primary purpose of the present study was to investigate the affects 
of flaws in unidirectional boron/aluminum composites, with the eventual 
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goal of predict ing the strength of such composites given a certain statistical 
distribution of internal flaws. These defects manifest themselves in two 
forms: a discontinuity in one or more boron fibers, or a localized void 

in the aluminum matrix. The loading condition of primary Interest is that 
of tension applied parallel to the fiber axes. With suitable modification, 
a so-called longitudinal model was analyzed with the micromechanics program 
in its original form to investigate the problems of modeling a flaw, generally 
a fibe.r discontinuity, and to evaluate the resulting localized stress concen- 
tration and the local plastic deformation it caused. The redistribution of 
the load to the broken fiber could also be characterized, but only up to the 
•Wnt at which a matrix element failed (crack Initiation). For further 
sti y of the load capability of the flawed composite, a crack propagation 
scheme is required. This capability permits a characterization of the 
energy required to isolate the defect in a "zone" of plastically deformed 
matrix material, or alternatively, the total energy capacity of the system 
at the point of catastrophic failure. 

The approach to crack initiation and propagation taken here is 
known as the "failed element" approximation as employed by Adams [4-6]. 

When an element in an area of high stress exhausts its strain energy capacity, 
It fails. From this, we assume that a "crack" has formed and has the 
dimensions of the failed element. This approximation has two implications, 
the most Important of which Is that a finite amount of material is removed 
from the system, which in an actual material is not the case. The other is 
that the crack is not likely to close up on itself in subsequent loading be- 
cause of its exaggerated width. These effects can be minimized to a practical 
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degree by making the finite element grid very fine and uniform in the 
area of anticipated crack initiation. 

It is not enough to simply delete an element from the finite element 
grid when it reaches its ultimate stress. The finite element method in- 
volves the maintenance of force equilibrium at every node point in the 
array, as discussed in Appendix A of the first-year report [1], This 
equilibrium must be maintained when an element fails or unloads. Thus, 
to represenc the unloading due to element failure, node point forces which 
are equal and opposite in sense to those equivalent to the state of stress 
within the element at its failure level must be applied at its node points. 

In addition, the failed element's material properties must be set to zero, 
so that the element makes no further contribution to the global stiffness 
matrix, and all of its computed values of stress and strain are set to zero. 
This insures that the element is completely unloaded and that no stresses 
will bo developed in it in subsequent load increments. 

In the present analysis, element failure can occur In one of two modes: 
when both the computed octrhedral shear stress and the plastic octahedral 
shear strain reach their maximum allowable values (maximum distortional 
energy criterion) , or when the hydrost tic tensile stress in an element 
exceeds the tensile ultimate strength of the material. This second failure 
criterion is also known as failure due to ultimate cleavage, and failure occurs 
whenever a tensile principal stress exceeds the. ultimate tensile strength. 
Complete details are given in Reference [1]. 



SUCTION 4 


AXTSYNMKTRIC ANALYSIS MKTIIOI) 

4.1. Purpose of the Axisymmetrie Formulatio n 

The development of an axlsyiranetrlc triangular finite element Cor the 
micromechanics analysis program was first proposed after the preliminary 
results of crack propagation studies using two-dimensional, generalized 
plane strain models had been examined. The difficulties and uncertainties 
of representing longitudinal loading of a square array of boron fibers 
embedded in an isotropic matrix led to the finite element models to be 
discussed in Section 6. As will also be discusse,.. in considerable detail 
in Section 6, the effect of fiber spacing on the state of stress and the 
pattern of crack growth around an Internal flaw is quite significant. It 
was thought that an axi symmetric formulation, allowing a fiber of circular 
cross section to be modeled, together with possible experimental correlation, 
might lead to the development of guidelines for the employment of the much 
more versatile generalized plane strain analysis. The parameter that is 
sought Is an '’effective" fiber spacing, that will accurately represent 
the relative affect of intact fibers surrounding a broken fiber, fibers 
that are, in general, not all equidistant from the flawed fiber. 

The axisymmetric finite element model provides a correct representa- 
tion of a single boron fiber JU? an annular sheath of aluminum matrix, as 
illustrated in Figure 1. The thickness of the annular sheath of aluminum 
can be changed quite readily, allowing the stress distribution around fiber 
flaws to be studied for various apparent fiber spacing conditions, with the 
circular cross section of the fiber being accounted for. It is proposed 
that a test specimen resembling this configuration be fabricated, and the 
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amount of displacement between the broken ends of the boron fiber at a 
specified load level be measured. This displacement could then be compared 
to the fiber displacement predicted by the mieromechanics analysis, and used 
to confirm its validity. In addition, by careful correlation of these 
results, including the state of stress predicted by the ax 1 symmetric formu- 
lation, to the strain and fiber displacement measurements of test specimens 
containing an arbitrary number of continuous fibers surrounding one dis- 
continuous fiber, an effective fiber spacing could be arrived at for use 
in the two-dimensional generalized plane strain analysis. The ultimate 
objective, of course, is to verify the limits of accuracy of the two- 
dimensional analysis, and so ovoid to the extent possible the much greater 
expense and complexity associated with using the three-dimensional formu- 
lation now available [7], 

Initially, the implementation of an ax i symmetric formulation was en- 
visioned to be a relatively simple task. However, due to the special 
nature of the existing micromechanics computer program and the requirement 
for maximum accuracy throughout the region being analyzed, the development 
of an uxisymmetrlc element proved to be considerably more complicated than 
was anticipated. Specifically, the manner in which boundary conditions and 
loading conditions are combined to aliow a unique solution technique 
necessitated the complete rewriting of these routines for use in the axi- 
symmetrio computer program. In addition, since the material lying on the 
axis of radial symmetry (the z-axis in Figure 1) must necessarily be included 
in the analysis, the familiar approximate form of the axisymmetric element 
stiffness matrix [12] was found to be unacceptably inaccurate. In the sub- 
sections that follow, the development of an exact axisymmetric triangular 
element is presented and the reasons for its special form are discussed. 



Fiber Discontinuity 



Figure 1, Region of Interest for the Axlsymmetric Analysis 
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4.2. The Axlaym metrlc Fi nite Element 
4.2.1. Basic AxlHvwmptrlc Re l atlonahipH 

The problem of analyzing a single fiber encased in an annular sheath 
of matrix material, and subjected to axial and radial tractions plus hygro- 
thermal gradients, falls into the class of problems known as torsionless 
axially symmetric states. These problems are generally defined relative 
to cylindrical coordinates (r,b,z), and can b« compared to the class of plane 
strain problems defined relative to Cartesian coordinates (x,y,z), as 
follows : 

o The in-plane stress components of the cylindrical system, 

o and o . correspond to o and a of the plane strain 
r z x y 

formulation, while o Q is the out-of-plane stress component, 

corresponding to a of the Cartesian system, 
z 

o The displacement components (u,v,w), corresponding to the 
(r,t\z) coordinates are such that (u,v) are independent 
of the polar angle, and the out-of-plane displacement component, 
v, vanishes. That is, u *■ f(r,z), w ■ f(r,z), and v ■ 0. 

When the strain-displacement relations and the stress-strain relations 
of the theory of elasticity are applied to this state of torsionless axial 
symmetry, it is found that “ r ^ * 0, and the stress components o^, 

0 rt , o and t are functions of coordinates (r,z) only. 

O’ z rz 

An axisymmetric finite element is in the form of a toroidal ring of 
constant cross section, as illustrated in Figure 2. The node points of 
such an element are in fact nodal circles, and the volume of such an 
element is dependent on both its cross-sectional area and the radii of these 
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FIGURE 2 . Solid Axisymmetric Triangular Element, 
Polar Coordinates 


nodal circles. In addition, nodal loads are a function of node point 
radii and the load per unit of circumference. The stress and strain vec- 
tors pertinent to Figure 14, as well as the strain-displacement relations, 
are shown below. Note that the out-of-plane strain at a point, e Q , is a 
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function only of the radial, displacement, u, and t Ho radial coordinate, r. 
Stresses ami et rains are related by a material properties matrix, [H], form- 
ing the following constitutive equation tor any given element: 

- flllji. > t (4) 

The material properties matrices are given in Appendix A far isotropic 
elastic and transversely Isotropic elastic materials, and for isotropic 
materials in the plastic range. In general, the forms of the material prop- 
erties matrices for axially symmetric conditions are identical to those 
found in the case of generalized plane strain, in that the coefficients 
relating the various corresponding stress and strain components are 
identical. However, most texts dealing with axially symmetric problems 
arrange the components of the stress, strain, and material properties 
tensors in an order different from that presented here. To have employed 
the more generally accepted sequence of arranging these tensor components 
would have made it necessary to rewrite all of those routines in the computer 
program in which stresses, strains, failure modes, and crack propagation 
are determined. To avoid this difficulty, all relationships are derived 
here with the stress and strain components arranged in the same order as 
found in the generalized plane strain relationships presented in Reference [1], 



4.2.2. Appro ximate? A xi symmetric. 


Ah is the* cast' with the generalized plane strain triangular element [1], 
a constant strain field Is assumed to exist within each element, which leads 
to the derivation of a shape matrix relating element strains to nodal dis- 
placements. While this derivation is much like that for a plane strain 
triangular element, the presence of the r-eoordlnate In the denominator of 
several of the terms leads to considerable difficulty in evaluating the 
coefficients of the. element stiffness matrix. The procedure is briefly 
outlined below, with Figure 3 provided as a visual aid in understanding the 
problem. 



FIGURE 3. Finite Element i in the r-z Plane 
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Bv ehoositiR a diBplacenrn.it' field in the simplest linear form (see 
Eq. A- l 2 and the related discussion in Appendix A of Reference {!]), wf 
arrive at the following relationship, 



where 


N 1 " 24 '“I + V + V 1 


N 2 ‘ 26 ln 2 + b 2 r + V 1 


N 3 ■ 24 |a 3 + V + c 3 al 


«>) 
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In which 


®1 " r 2 Z 3 ” r 3 Z 2 


r 3 Z l ‘ r l Z 3 


a 3 ’ r l Z 2 ‘ *2*1 


b l " Z 2 “ z 3 


b 2 « *3 > \ 


b 3 “ Z 1 z 2 


( 8 ) 


r 3 ” r 2 


C 2 " r l " r 3 


r 2 " r l 


the 4x6 rectangular matrix in Eq. (7) is the "shape matrix," [B] , for the 
axisymmetric triangular element, and can be used to form the element 
stiffness matrix for Individual elements. Note that the equation for the 
shear strain, y, , has coefficients with r terms in the denominator. When 

I Z 

a node (or nodes) of any element lies on the axis of rotation, 1 . e . , r « 0. 
singularities in the shape matrix result. In addition, when one considers 
the operation of forming the element stiffness matrix, 

[k] -/ [B| T [H] [B]d (Vol ) (9) 

* Vol 

and the form of the shape matrix, [B], in Eq. (7), it is obvious that a 
term-by-term integration involving the r and z coordinates of the node 
points is necessary. The r terms that appear in many of these expressions 
result in some rather tedious calculations, and lead to logarithmic terms 
which can also result in singularities in the stiffness matrix. In the 
case of a plane strain or generalized plane strain clement, the volume 
integral required by Eq. (9) is simply equal to the cross-sectional 
area of the element times its thickness, as described, for example, in 
References [11,12], To avoid the difficulties presented by the r terms 
in the denominator of the shape matrix terms, an average shape matrix, [B] 
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can bo formed, using the coordinates of the centroldal point of the triangle 


as the element's coordinates. 

and an element stiffness matrix can be calcu- 

lated directly, i.e., 


[kJ 1 - 27r[B) T [H](B]FA 

(10) 

where 


r " (r t + r 2 + r 3 )/3 

(11) 

and 


z - (z L + z 2 + z 3 )/3 

(12) 


Solutions using the approximate element stiffness matrix of Eq. (10) have been 
found to be quite acceptable as long as the planar dimensions of the indivi- 
dual elements are small compared to their radial coordinates, say on the 
order of 10 to 1. In particular, hollow cylindrical bodies can be very ade- 
quately analyzed using the element stiffness formulation of Eq. (10). How- 
ever, when the cylindrical body Is solid, or possesses a very thick wall, 
large numerical errors are encountered when the approximate element stiffness 
matrix is used. In the present analysis, the axis of rotation is necessarily 
part of the region of investigation, and being composed of the very stiff 
boron fiber, carries a significant portion of any applied axial loads. Trial 
analyses using the approximate ccntroidal formulation of Eq. (10) indicated 
that very large numerical errors, on the order of 15 to 20 percent, were 
present In stress components that were normal to the direction of :he applied 
axial load. After examining these numerical errors for several loading 

situc-ions, and considering the need for the precision required if the 

4 

objectives of the axisymmetrlc analysis were to be met, it was decided to 
formulate an exact axisymmetric element stiffness matrix, i.e., to perform 
the integrations indicated in Eq. (9). 
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4.2.3. Exact Axisymmetric Element Formulation 

In this section, a fairly thorough description of the derivation of 
the integral form of the axisymmetric triangular element is presented. In 
order to maintain the flow of the derivation, some of the more involved 
mathematical procedures have been placed in appendices, while the results 
of the procedures are used directly. 

4 . 2 . 3 . 1 . Strain - Displacement Relationships 

As is the cose for planar, constant strain finite elements, a simple 
linear displacement field is chosen for the axisymmetric triangular element. 
This classifies it as a "simplex" element, but unlike the planar, or unit 
thickness, elements, it is derived in terms of the element's generalized 
coordinates, {£}. Planar elements, with their simpler strain-displacement 
relationships, allow direct evaluation of the shape function coefficients in 
terms of the finite element model's global coordinates. The displacement 
field relating the displacements of a point in the region of analysis to its 
generalized coordinates for a polar problem can be expressed as. 



This relationship can be shown to meet convergence criteria for conforming 
finite element displacement fields [10]. Substituting the nodal coordinates 
of an element such as that shown in Figure 3 into Eq. (13) yields 
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(14) 


or 

- [T) i (U 1 (15) 

where {6}^ is the element nodal displacement matrix, {f,}^ Is the element 
generalized coordinates matrix, and [T]^ is the element transformation matrix 
that relates the two. By examining the transformation matrix in Eq. (14), 
it can be seen that only translations are involved. In other words, the 
element is not derived with respect to a "natural* 1 coordinate system, as 
in the case of a beam or isoparametric element. However, due to the need 
to evaluate specific integral coefficients, the terms f,, through £, cannot 
he evaluated directly in terms of the global coordinates, as is the case for 
a plane strain simplex triangle. Solving Eq. (15) for ^ wo obtain an 
equation relating the generalized coordinates to the nodal displacements, 


- ' T k lls) i 


( 16 ) 


where the inverse of [T]j can be shown to be [3] 
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tn which A equals the triangular cross-sectional area of the element, and 
the coefficients a^, b^ and c^ are as defined In Eq. (7). 

Substitution of the assumed displacement field relationships, Eq, (14), 
into the definition of the strain components, Eq. (2), leads to 
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(17) 


or In matrix notation. 


Ml . [Cl^.1, 


(18) 


wliero [C] ± is the shape matrix in terms of the element generalized 
coordinates. For global coordinates and displacements, we combine Eq. (16) 
and Eq. (17) to obtain 


tr h ■ tC] i tT,- l UJ 1 

which is the desired strain-displacement relationship. 


( 19 ) 
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A.2.3.2. Element Stiffness Matrices 

For the analysis of highly stressed, fiber-reinforced, metal matrix 
composites, three distinct element stiffness matrices are required; o ie to 
model an isotropic matrix material which is loaded below its elastic limit, 
another to model the fibers, which may be transversely isotropic, and a 
third to describe the behavior of the isotropic matrix when it is loaded 
into the plastic region. These three classes of material response are 
discussed in detail in Appendix A of Reference [1], These same relation- 
ships are used to derive the required axisymmetric element stiffness matrices. 

In general, a stiffness matrix relates nodal forces to nodal displace- 
ments. With the forces as known quantities, this allows the nodul displace- 
ments throughout the region of analysis tp be solved for, and stress 
components can then be "backed out" of these displacements. The integra- 
tion of the product of the shape matrix, [C]^ and the stress-strain matrix, 
[H]^, over the volume of an element yields its element stiffness matrix, i.e., 

[«,./ [B]J[H] 1 [B] 1 d(Vol) (20) 

Vol 

or, for an element stiffness expressed in generalized coordinates, 

[k ] 1 -J [C]^[H] ± tC] 1 d(Vol) (21) 

Vol 

For an elastic, isotropic material, after substituting Eq. (17) and Eq. (A-14) 
of Reference [1] into Eq. (20) and performing the indicated matrix multi- 
plications we have, after integrating with respect to 6 , 
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Integrating Eq. (22) with respect to r, 
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where E ■ Modulus of elasticity 
V ■ Poisson's ratio 
and the integrals are represented as 
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*i ’ \j z IiTdr - x « " J c l z ^ drdz 

*2 ■ IJ, drdl '* ’ U^ rdS 

h'U Zdtdl *« ‘ IJ^ rdl 

r z 


(24) 


The first three integrals are easily evaluated, and are defined in several 
finite element analysis texts, e.g., Reference [11]. They are 


(rj+r 2 +r 3 )[r 1 (z 2 -z ;} )+r 2 (z 3 -z 1 )+r 3 (z :i --z 2 )] A(r 1 +r 2 +r ;} ) 

X 1 ” ~~ 6 " 3 

(r 1 (z 2 -z 3 )+r,(z 3 -z 1 )r 3 (z 1 -z 0 ) ] 

I 2 = ~~2 ' A 

(Zj+z^zp [r 1 (z 2 -z 3 )-fr 2 (z 3 -z 1 )+r 3 (z 1 -z 2 ) ] A(z 1 +z 2 +z 3 ) 

T m — — 1 ' ■ - ' 1 - “""" 

3 6 3 


(25) 


The integrands of I ^ , I,. and 1^ contain a (^) term, and are considerably 
more difficult to evaluate. 

In addition, when r^ , r 9 or r 3 lie on the z-axis, the integrand be- 
comes singular and special procedures must be employed to evaluate these 
terms. The procedure for evaluating these Integrals is presented in 
Appendix A, and the manner in which singularities in the integrands are 
dealt with is described in Appendix B. For the most general rase, however, 
the expressions for these integrals are presented below. For 



each integral, the symbol L, is used [11] to indicate that a cyclic permu- 
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tat Ion of the indices is necessary to obtain the full expression, l.e.. 
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In a similar manner. 
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The element stiffness matrix for the transversely Isotropic (fiber) 
case is obtained in exactly the same manner as for the isotropic elastic 
ease, except that the material properties matrix presented in Eq. (A-16) 
of Reference [1] is substituted into Eq. (20). This leads to the following 
transversely isotropic, generalized, element stiffness matrix, 
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where Q - (l+v)(l-\>- - ^ - V ; ) 

L 



F - v'(l+v) 


E' - Elastic modulus in the direction perpendicular to the plane 
of isotropy 

v' * Poisson's ratio representing a strain in the plane of 

isotropy due to a normal stress in the direction perpendicu- 
lar to that plane. 

For the element stiffness matrix of an isotropic material in the plastic 
range, the material properties matrix presented as Eq. (A-32) in Reference 
[1] is used in Eq. (20). After performing the required multiplications and 
integrations, we have, in terms of generalized coordinat ^s, the isotropic, 
plastic element stiffness matrix. This expression is shown in Eq. (30). 

The terms are components of the deviatoric stress tensor, t q is the 



3t [1+C1+v)1 


23 


octahedral shear stress and 2Hj, Is the tangent modulus of the octahedral 
shear stress-octahedral plastic shear strain curve, all of which are discussed 
in detail in Appendix A of Reference [1]. 

Each of the element stiffness matrices given by Eqs. (23), (29), and 
(30) must be expressed in terms of the global coordinates of the finite 
element model before they can be used to assemble the global stiffness 
matrix. This is accomplished by evaluating the inverse of the transformation 
matrix for each element, as defined in Eq. (16), and using the following 
relationship : 

[k] ± - (rT]” 1 )^[k) i [T]^ 1 (31) 

The global stiffness matrix, [K] , is assembled, element by element, by 
a subroutine which also imposes the boundary conditions required for the 
specialized loading technique (see Appendix A-4 of Reference [1]). In this 
subroutine, each element in the model is examined to determine whether it is 
fiber or matrix, elastic or plastic, and the appropriate element stiffness . 
subroutine is called. In these element stiffness subroutines, the strain- 
displacement matrix and the stress-back substitution matrix for each element 
is evaluated and stored on a peripheral device. These are required to oh;ain 
element strains and stresses from the nodal displacements that the finite 
element solution provides. In general, 

(o} i - [H] i [C] 1 [T}" 1 1 {6} 1 - [H] 1 [A] i (6} 1 - [B] ± {6} i (32) 

where [A]^ is the strain, -di. ’nlacement matrix and [B]^ is the stress-back 
substitution matrix. 

These matrices have to be evaluated for each of the three material 
conditions. The most general form of the [A] ^ and (B]^ matrices are shown 
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where T, !■’, and <) are as defined in Eq. (29). The stress-back substitution 
matrix for a plastic isotropic material is as shown in Eq. (35) on the next 
po^e, with the terms A, \’, B, and as defined in Eq. (30). It is 
interesting to compare the elastic and plastic cases, especially in the 
terms that are zero for the elastic [Bj matrix and negative for the plastic 
matrix. As is also the case for the plastic element stiffness matrix 
versus the elastic form, the negative sense of the additional terms is 


a reflection of the reduced module© of most elastic-plastic materials when 
in the plastic range. 

The strain-displacement matrix, being a function of element geometry 
only, is the same for all three material conditions, i.e., 
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(36) 


It is important to note that for a given set of nodal displacements, 
Eqs. (33), (34), and (35) describe the variation of stresses within the 
element as a function of the r and z coordinates. In other words, using 
these relationships, the exact state of stress at any point in the plane of 
the element c a. he obtained. In this study, the centroid location of each . 
element has been nrogrnmmed in as the element point of interest. 

This concludes the description of the exact, triangular, constant 
strain axisymmetric element. It is important to realize that those elements 
that have a node point or a side coincident with the axis of rotation, 
sometimes referred to as "core" elements [12], require special treatment. 
Elements having one side parallel to the z-axis also lead to singularities. 
These special cases, described in detail in Appendix B, make the computer 
implementation of this formulation particularly complicated, as the 
strain-displacement and back substitution matrices must also be modified. 
Comparison of the exact formulation of the element to the approximate form 
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indicates that the additional effort necessary to develop and Implement 
the Integral form Is justified. For axial loads applied to a model of 
"core" elements, the error in a stresses is on the order of 14 percent for 
the approximate formulation. The error in o Q stresses is on the order of 
18 percent. For radial or combined loads, the error observed in the 
approximate element formulation is even larger. 



SECTION 5 


MATERIAL PROPERTIES 


In modeling the boron/ aluminum composite, the boron fibers have been 
treated as brittle, linearly elastic materials with isotropic strength 
and stiffness properties. The aluminum matrix has also been considered to 
be isotropic, but is modeled as an elastoplastic material. To accomplish 
this, the actual stress-strain curve of the aluminum alloy selected is 
input to the analysis by curve fitting via a Richard-Blacklock two-parameter 
equation [13], as discussed in Appendix A-5 of Reference [1], Thus, at any 
load level the tangent modulus for any given element can be computed. This 
makes possible an accurate representation of the plastic deformation of 
the matrix. 

Although the nonlinear material properties of any matrix material, 
e.g., another aluminum alloy, can readily be incorporated in the analyses, 
a 6061-T6 aluminum alloy at 75°F was used in obtaining the present results. 
The material properties shown in Table 1 were obtained from Reference [14]; 
the full range stress-strain curve for determining the curve fit parameters 
used is shown in Figure 4. 

TABLE 1 

Aluminum Matrix Material Properties - 6061-6 Alloy [14] 

£ 

Young’s Modulus E = 10.0 x 10 psi 

Poisson’s Ratio v =0.33 

Tensile Yield Strength F^ * 36000 psi 

Tensile Ultimate Strength F tu = 45000 psi 

Coefficient of Thermal Expansion a = 13.0 x 10~ / °F 
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The boron fiber properties indicated in Table 2 were obtained from 
Reference [15]. 


TABLE 2 

Boron Fiber Material Properties [15] 


Young's Modulus 
Poisson's Ratio 
Tensile Ultimate Strength 
Ultimate Strain 

Coefficient of Thermal Expansion 


E ■ 60.5 x 10^ psi 

v - 0.130 

ptu « pty . 500,000 psi 

ptu - F tu - 8.264 x 10" 3 in. /in. 
E 

a - 9.0 x 10-6/°F 



Figure 4. Typical Full Range Stress-Strain Curve for 6061-T6 Aluminum 
Alloy at Room Temperature [14] 


SECTION 6 


NUMERICAL RESULTS 

In this section, the finite element analysis methods discussed in 
Sections 3 and 4 are applied to a variety of finite element models. 

In Sections 6.1 and 6.2, the development of the various finite element 
models is discussed, for use with the generalized plane strain analysis 
(described in Section 3) and the axisymmetric analysis (described in Section 
4) , respectively. 

In Sections 6.3 and 6.4, numerical results are presented for axial 
loading of the longitudinal section models, while axial loading of the 
transverse section models is presented in Section 6.5. Transverse loading 
of the transverse section model is discussed in Section 6.6. 

In Section 6.3, the results of the axial loading of models representing 
a condition of 33 percent discontinuous fibers are discussed. At the time 
these results were generated, early in the present second-year study, the 
basic computer program, developed during the first-year study [1], was still 
in a somewhat unrefined condition; nearly all the reduction of output data 
had to be done by hand. With the conversion of the computer program to 
the larger, faster Control Data Corporation Cyber 730/760 computer system, 
installed at the University of Wyoming in early 1980, a post-processing 
package was added which is capable of drawing the material interfaces of the 
finite element model, the outline of any crack that might be present, and 
a variety of stress and strain contours, as specified by the user. The 
results of the analysis of the considerably larger and more complex models 
then possible are discussed in Section 6.4. 
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Sections 6.5 and 6.6 deal with the mechanical loading of transverse 
section finite element models. As a result of these studies, transverse 
section models have been found to be Incapable of adequately representing 
stress concentrations due to material defects when leaded in the direction 
of the fiber axes. However, they are especially useful in studying the 
effects of transverse mechanical loads and hygrothermal loads. 

Preliminary results of the newly developed axisymmetric analysis are 
presented in Section 6.7. These results demonstrate the capability of the 
analysis; more investigation remains to be done. 

6.1. Generalized Plane Strain Analysis Models 

6.1.1. Development of the Broken Fiber, Longitudinal Section Models 

There are two primary reasons for the development of a longitudinal 
section model. One is to permit study of localized stress concentrations, 
the resulting elastoplastic behavior of the aluminum matrix, and subsequent 
crack propagation in the area of fiber flaws. Another is to characterize 
the load carrying capability of a flawed fiber as a function of distance 
from the location of the fiber flaw. 

These two considerations lead to the most important aspects of designing 
the longitudinal section models, i.e., geometry, finite element grid 
resolution, boundary conditions in the vicinity of a flaw, and spacing 
of the boron fibers in the model. The problem of fiber spacing will be 
discussed first. 

A typical cross section of a unidirectional, square array, boron/ 
aluminum composite as shown in Figure 5, the section being perpendicular 
to the fiber axes. A longitudinal section finite element model attempts 
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FIGURE 5. Cross Section of a Square Array of Fibers, 
55 Percent Fiber by Volume 


to represent the composite in a plane oriented perpendicular to this 
section. A longitudinal model of a section parallel to the x or y axes, 
through the centers of the fibers, would be representative of the minimum 
distance between fibers. A section cut at 45" to the x-axis and through 
the F iber centers would depict a maximum fiber spacing situation. When one 
of these fibers in broken, the load it carries decays to zero at the broken 
surface, assuming that the boron-aluminum interface remains Intact. At the 
flaw site, the fibers adjacent to the broken fiber, and to some extent the 
surrounding aluminum matrix, must absorb the load that the broken fiber would 
have otherwise carried. The aluminum transfers this excess load back into 
the broken fiber via a shear mechanism so that at some distance from the 
fiber break, that fiber is again fully effective in carrying load. It is 
logical to presume that the amount of aluminum between the boron fibers will 
have an effect on this load transfer mechanism. To characterize the effects 
of variation in fiber spacing, two longitudinal models were studied, one 
representing a 90° section cut of the transverse cross section, and another 
representing a 4^° section cut. A 45° section model is shown in Figure 6, 
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the element numbers being indicated. The region In the lower left corner, 
with unnumbered elements, Is the region in the vicinity of the broken fiber 
end. This local region is shown in the expended view in Figure 7, the 
clement numbers being given here. A 90* section model is shown in Figure 8. 
Note that the fiber diameter dimensions have been normalised to unity. In 
Figure 8 the effect of the 90* section cut in diminishing the amount of 
local aluminum matrix la shown quite clearly. The else and aspect ratios of 
the fiber elements are exactly the same as those of Figure 6, but the 
aluminum elements of the 90* section model are so compressed that the 
element numbers, which ore identical to those of Figure 6, have been 
eliminated for clarity. 

The second problem that must be s^ved in the finite element modeling 
of a broken fiber in a composite is the geometry in the area of the fiber 
discontinuity-matrix interface. The efforts of the present investigators 
to resolve this problem have been evolutionary in nature and many models 
were developed and discarded in the process. It will be noted, for example, 
that the models of Figures 6 through 8 are slightly different in the region 
c the broken fiber end than the modela presented in Figures 3 6 through 
3.8 in Reference [1]. 

As the longitudinal finite element models evolved, the problem of 
computer capacity came to be a limiting factor. The models shown in 
Figures 6 through 8 yielded generally favorable results, but represented 
the limits of the capability of the Xerox Sigma 7 computer then being used. 
After study of the numerical results obtained using these models, it was 
decided that they were limited in three important areas. The first was 
the modeling of the region in which crack growth is expected to occur. 
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FIGURE 7. 


Finite Element Representation of 
Broken Fiber End in Figure 6. 


the Local Region Around the 
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Ideally, this region should contain an extensive area of uniform elements 
of approximately the same size as the elements forming the initial crack. 

This is due to the fact that hen a significantly larger element is en- 
countered by the growing crack, propagation ceases until this larger element 
and its neighbors have been strained to their ultimate value. When such a 
large element falls, the released energy to be redistributed among the 
surrounding elements is considerable, and the result is that many more 
elements fail in the process. This also dictates the formation of a fairly 
large cavity within the model, which may affect the pattern of subsequent 
crack propagation in an unrealistic munner. Tills was particularly noticeable 
in the case of the 90° section model. A second limitation of a small model 
is that it represents a situation In which fully one-third of the fibers 
are broken. In other words, the effects of a broken fiber on the loading 
of more remote Intact fibers cannot be studied; fibers which may have a 
considerable affect on the pattern and extent of crack growth. Finally, 
the limited axial length of the small models of Figures 6 and 8 was thought 
to be an unfavorable influence in terms of end O-ffects. As will be discussed, 
this limitation prevented the loading of the material being modeled to 
its full capacity due to the arrival of the crack front at the right boundary 
of the model. 

With the acquisition of the much larger and faster CDC Cyber 760 
computing facility, large: and more complex models could subsequently be 
studied. The two larger finite element models created for this purpose 
are presented in Figures 9 and 10, representing 45° section and 90° section 
longitudinal models, respectively. In them, 12.5 percent of the fibers 
are discontinuous. It will be noted that, in addition to an extensive 
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FIGURE 10. Refined 90° Section Longitudinal Model wit 12.5 Percent Discontinuous Fibers, 
8.2 Fiber Diameter' in Length. 



region of uniform matrix element!? adjacent to the filter discontinuity, 
many rows of elements are retained in the region between this fiber and its 
neighbor, along the full length of the model. 

6.1.2. D evelop ment of the Transverse Se ct ion. Models 

Finite element modeling of i transverse section of a unidirectional 
boron/aluminum composite Is fairly straightforward [4—6 ,8,9] . A typical 
transverse section model is shown in Figure 11. This model will be used 
here to demonstrate crack propagation In an unflawed composite subjected 
to a transverse normal loading. However, the need to study the influence 
of a reduced load capacity in one fiber on Its neighboring fibers requires 
that a minimum section model such as that shown in Figure 12 be employed. 

This model represents the first quadrant of a repeating square array of four 
fibers. If the fiber centered at the origin is assumed to be a flawed 
fiber, it will be surrounded by eight other (unflawed) fibers in the array. 

A model of this type can easily lead to a great number of finite elements, 
and attempting to increase the resolution of the grid at selected locations 
often results in a very large bandwidth of the overall stiffness matrix 
for the finite element model. 

6.2. Axlsymmet r ic Ana lv sis Models 

The finite element model used for the preliminary studies using the 
axisymmetric element formulation is shown in Figure 13. In this model, the 
horizontal axis is the r-axis, while the vertical axis is defined as the 
z-axis, or axis of rotation. The fiber elements are located along the z-axis 
and extend radially outward for three "bays" of elements. The fiber 
discontinuity is modeled by freeing the first four node points at the lower 
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FIGURE 12. Transverse Section Finite Element Model of a Repeating Square 
Array of Four Fibers, 55 Percent Fiber by Volume. 
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Matrix 


r ,u 


FIGURE 13. Axisymmetr 1c Finite Element Model Representation of the 

Maximum Fiber Spacing in a 55 Percent Fiber Volume, Square 
Array, Unidirectional Composite ( r f/ r m * 0.424). 

Model Length is 5.5 Fiber Diameters. 
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loft corner of the model, extending radially outward. The philosophy in 
creating this model Is similar to that followed in the specificiation of 
the generalized plane strain models of Section 6.1, i.o., maximum element 
resolution in areas of high stress gradients. Minimum stiffness matrix 
bandwidth is achieved by maintaining continuous node point numbering from 
one end of the model to the other in at least one coordinate direction. 

This can lead to a few more elements than are required in some cases, but 
it has been found that the savings in computer core space due to a reduced 
bandwidth far outweigh the cost of a few superfluous elements. In 
addition, the dimensions and matrix thickness of this model are very easily 
changed, making it particularly useful for parametric studies of the effects 
of fiber volume on the response of this particular configuration. 

t> . 3 . Axia 1 Loading of Longitudinal Models With 33 Percent Discontinuous 
Fibers 

6.3.1. Crack Initia tion and Propagation in the 45 ° Secti on Longitudinal 
Mod ex 

The 45° section longitudinal model was loaded axially with one boron 
fiber treated as discontinuous, the broken ends being in contact when 
loading was Initiated. Crack initiation occurred with the failure in 
octahedral shear stress of Element No. 2 (see Figure 7) at an average 
applied stress of 64.1 ksi. The release of energy caused by the failure 
of Element No. 2 resulted in the failure of seven additional elements. 

Nos. 1, 3, 4, 13, 14, 15, and 26. With the stiffness capacity of these 
elements deleted from the analysis and their strain energy redistributed to 
the remaining model, 13 more elements failed. This process of crack growth 
at constant stress continued in a regular pattern until the "crack" had 



progressed to the point shown in Figure 14. The elements that have failed 
are blacked out. 

After crack growth had ceased, monotonlc loading of the composite was 
continued to an average applied stress level of 107.2 ksl without further 
element failures. As loading progressed, the aluminum matrix elements 
adjacent to the fiber nearest the discontinuous fiber experienced increasing 
amounts of plastic straining. The shaded elements in Figure 14 are those 
in the plastic stress range at an average applied stress level of 107.2 ksi. 
An examination of the in-plane components of stress for these elements 
revealed that in the elements nearest the crack tip, shear stress was of 
the greatest magnitude. In the plastically strained elements farthest 
from the crack tip, tensile stress, parallel to the fiber axes, was again 
the major stress component, although the shear stress level was still high. 

As the crack formed and grew, the load level in the broken fiber 
decreased relative to that of the intact fibers, as expected. The pattern 
of crack growth exhibited by this analysis, and the manner in which each 
element was deformed and failed primarily by shear stress, is very similar 
to experimental results obtained by Awerbuch and Hahn [16]. In their 
study, center-notched tension specimens of unidirectional boron/aluminum 
were tested. Microscopic examination of the failed test specimens revealed 
crack growth in the aluminum matrix adjacent to the last cut boron fiber 
on either side of the notch. These cracks appeared to propagate parallel 
to the fiber ,es, and were accompanied by long zones of plastic shear 
deformation in the matrix, also running parallel to the fiber axes. 

At an average applied stress of 107.2 ksi, Element No. 136, located 
between the two continuous fibers (see Figure 6), failed in hydrostatic 



ORIGINAL PAGE 13 
OF POOR QUALITY 


uoj 


xmjpw 
uinujuin jv 


FICURF 14. Extent of Crack Growth and Plastic Deformation at an Average Applied Stress of 

107.2 ksi, 45° Section Longitudinal Model. Crack Initiation Occurred at 64.1 ksi 







tension. Subsequent load redistributions resulted in total fuilure of the 
aluminum matrix. At this load level, the fiber adjacent to the discontinuous 
fiber was carrying a great deal more load than it would have were there no 
broken fiber. The aluminum matrix, loaded primarily due to strain compat- 
ibility with the boron fibers and by Poisson effects, is predicted to rupture 
when the. analysis is confined to the boundaries of these smaller finite 
element models. For a longer finite element model with a higher percentage 
of continuous fibers, higher average applied stress levels can be sustained 
without failure, as will be shown loter, and further growth of the crack will 
occur. 

The composite stress-strain response is plotted in Figure 15, as a 
measure of the strain energy capacity of this 45° section boron/aluminum 
composite in which one fiber out of three is discontinuous. It will be 
noted that the rate of energy absorption with increasing stress after 
crack initiation is considerably greater than that exhibited up to the point 
of initial failure. That is, the slope of the stress-strain curve is less. 
From this plot it is obvious that plastic deformation and crack growth are 
Important considerations in the evaluation of the effects of flaws in 
composite materials. 

The large amount of straining of the total model that takes place 
during crack formation and propagation will also be noted in Figure 15. This 
further illustrates the effect and extent of the crack growth illustrated 
by Figure 14. Finally, the loss of some of the broken fiber's effectiveness 
in carrying the applied load is clearly illustrated in Figure 15 by the 
significant reduction in the composite modulus after crack formation. This 
change in modulus is particularly dramatic in this case because of the high 
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FIGURE 15. Plot of Applied Stress versus Composite Strain for the 45° Section Longitudinal 
Model, 33,3 Percent Discontinuous Fibers 
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percentage (33 percent) of discontinuous fibers, the extent of crack propa- 
gation, and the relav short axial dimension of the model being investi- 

gated. 

6.3.2. Crack Init iation a nd Propagation in the 90 ° Section Longitudinal 
Model 

When the 90° section longitudinal model was loaded in a direction 
parallel to the fiber axes, with one fiber broken, a large stress concen- 
tration occurred at the crack tip, as was the case with the 45° section 
longitudinal model. However, the composite axial stiffness of the 90° 
section model is considerably greater than that of the 45° section model, 
due to the larger fiber volume fraction of the 90° model. 

The closer proximity of an intact fiber with th broken one in Che 90° 
model results in a greater shear stress gradient at the end ol the broken 
fiber than was seen in the 45° model. As a result, crack initiation occurs 
at an average applied stresb level of 37.1 ksi, due to the failure of 
Element No. 2, as defined in Figure 7. Subsequent load redistributions re- 
sulted in a series of a single element failures, until finallv t' ’lure 
of Element Nv/. 42 triggered the failure of Element Nos. 53, 54, and 55, 
with Element Nos. 65, 66, and 67 failing after that. At this point, crack 
propagation ceased, with plastic deformation around the crack tip and along 
the fiber progressing as the load level was increased to 84.0 ksi. Figure 
16 illustrates the pattern of crack growth and plastic deformation at this 
level of applied stress. 

Tn comparing Figure 16 with Figure 14, the difference in the shape 
and extent of the crack is as obvious as the difference in load level. 

The higher shear stress gradient brought about by the closer fiber spacing 




FIGURE 16. Extent of Crack Growth and Plastic Deformation at an Average Applied Stress of a 
ksi , Q 0° Section Longitudinal Model. Crack Initiation Occurred at c = 37.1 ksi 
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causes earlier crack initiation, but the crack appears to be promptly con- 
tained by a large zone of plastically deformed matrix. A strong similarity 
between the two is shown by the puttern oi plastic deformation and in the 
fact that it is due primarily to shear, particularly In the vicinity of the 
crack tip. Again, the plastic deformation appears to progress down the 
boundary of the unflawed flbe!r adjacent to the broken fiber, which is con- 
sistent witli the pattern observed in the 45° model, and with the experiments 
of Awerbueh and Hahn [16]. 

At an average stress level of 84.1 ksi. Element No. 78 failed (see 
Figure 7), and subsequent element failures resulted in the crack pattern 
shown in Figure 17. At this point, the "crack" was over 5 fiber diameters 
long, extending to the opposite boundary of the finite element model. The 
broken fiber Is carrying very little load under this condition, as it is 
now only connected to the remainder of the model by a single node point at 
the right boundary. With continued loading, the last of the matrix elements 
adjacent to the broken fiber failed at 98.0 ksi, and their release of 
energy triggered the rupturing of Element No. 134 (see Figure 6) by hydro- 
static tension. 

The energy absorption capacity of Che 1 lawed 90° section model of a 
boron/aluminum composite is indicated by the stress-strain plot of Figure 
18. Very little pure straining takes place, in the initial stage of crack 
growth, as depicted by Figure 16. Both Figures 15 and 18 are plotted 
to the same scales of stress and strain so that the composite axial stiff- 
nesses of the 45° and 90° section models may be compared visually. As 
was the case with the 45° section model, the 90° section model exhibits 
a reduction of stiffness after crack initiation and propagation, due to a 




90° Section Longitudinal Model. 




AVERAGE APPLIED STRESS, c (ksi) 
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change Is l.ir I t a Hs pronounced in tlu> 90° section model, .is would he expected 
with its nuu'h smaller initial crack formulation. 

o.i. Axial Lomling of Long itudlna 1 Mode is with 12.5 Per cent Dis continuo us 
F i hers 

Crack In it l.it ion and Propagat ion I n the 45° S ection Longitudina l 

Mode 1 

The 45° section longitudinal model was loaded axially with one boron 
tilier again treated as discontinuous. Plastic deformation around the stress 
concent rat Ion caused bv the fiber discontinuity was first observed at an 
average applied stress level of 18.0 ksi. Loading was continued until an 
initial failure occurred at an applied stress level of 74.1 ksi, which is 
considerably higher than the 64.1 ksi level at which first failure occurred 
in the 44° mode! studied in Section 6.3.1. In addition, only one element 
tailed, the crack tip being temporarily "blunted" by this failure. The 
state of stress in the aluminum matrix Just prior to the initial failure 
is clearlv illustrated in the contour plots in Figures 19 through 22. 

Although anv or all of the various stress and strain components can be 
plotted by the computer program, contours of constant octahedral shear 
stress, octahedral shear strain, maximum principal stress, and in-plane 
suear stress have been selected here as being the most useful in studying 
the axial loading of unidirectional composites containing defects. Note 
that the octahedral shear stresses in Figure 19 have been normalized with 
respect to the octahedral shear yield strength of the 6061-T6 aluminum 
allo\. In this wav, the region of plastic deformation can be readily 
discerned, as any contour value equal to one defines a plastic zone boundary. 




Section Longitudinal Model. 
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FIGURE 20. Contours of Constant Octahedral Shear Strain (10 in/in). Average Applied 
Stress,^ * 75.1 ksi, 45° Section Longitudinal Model. 
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It will be* noted that at an average applied stress level of 75.1 ksl, the 
plastic zone Is on the order of one fiber diameter In length. The state 
of atreaa in the flbera Juat prior to crack Initiation la alao of great 
Interval ; thia la llluatrated quite effectively in Figure 23. In thia 
figure, the atreaa level of each of the flbera in the model la plotted aa 
a function of dlatance from the alte of the fiber dlacontinulty. Aa 
expected, the broken fiber plcka up load fairly quickly via ahearlng atreaaea 
in the matrix, attaining £9 percent of the load level of the moat remote 
flbera at 8.2 fiber dlametera (end of model) from the plane cf the dia- 
contlnuity. Of apodal intereat la the affect the broken fiber has on lta 
nenreat neighbor, and to some extent, the next fiber alao. Theae flbera 
must make up for the loat load carrying capacity of the compoaite in the 
region of the fiber dlacontinulty, with the adjacent fiber carrying 10 
percent more load than the remote fibers. It will alao be noted that the 
remote fibers gradually Increase in loading as a function of distance from 
the site of the fiber break, while the two closer fibers tend to unload. 

At an average applied j stress of 76.0 ksi, crack growth involving 
the failure of 82 elements took place. This process required 16 redistribu- 
tion steps, i.e., the failure of one group of elements would trigger the 
failure of additional elements, causing the crack to continue to propagate 
at the constant load level. The state of stress in the matrix Just after 
this period of crack growth is shown in Figures 24 through 27. Note that 
the size and shape of the crack is plotted as well as the contour values. 

Note also the reduction in the size of the plastic region, as ^hown in 
Figure 24, and the increased stresses in the matrix between the other 
fibers, especially the shear stress. The fiber loading plot of Figure 28 
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FIGURE 26. Contours of Constant Maximum Principal Stress (ksi). Average Applied Stress 
a * 76.0 ksi, 45* Section Longitudinal Model. 
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corrcs ponds to the state of stress and deformatior of Figures 24 through 27, 
and the most change caused by the element failures appears to be the redueed 
loading of the broken fiber. While It Is still loaded to almost the same 
percent at the right boundary of the model, the slope of its loading curve 
is somewhat less near the break site. 

The next set of figures, Figures 29 through 13, are pi >ts corresponding 

to an average applied stress of 97.4 ksl. At this load level, further crack 
growth was Immanent. It will be noted that the region of plastic deforma- 
tion has enlarged considerably, and that the loading of the broken fiber 
at 8.2 fiber diameters from the plane of discontinuity is now only 86 percent 
of the load carried by the two most remote fibers, Fibers 4 and 5. An increase 
In the load resulted in only one additional element failure, however, and 
cr ick propagation did not resume* until an average stress level of 183 ksi 
had been attained. 

At an applied composite stress of 183 ksi. crack propagation involving 
169 elements and 11 stages of growth at constant stress occurred. The 
results of this growth are illustrated in Figures 34 through 38. What is 
most significant about this stage of the loading of this model is that tin* 
extensive growth of the primary crack has resulted in such a high, localized 
stress in the adjacent fiber that two large matrix elements between it and 
the next intact fiber have ruptured. As Figure 38 illustrates, the adjacent 
fiber (Fit :r 2) is now carrying about 14 percent more stress than Fibers 4 
and 5, and appears destined to fall well before Fibers 3, 4, and 5. 

The 43° section longitudinal model continued to absorb loading up to an 
average stress level of 208.5 ksi with no further failures. The analysis 
was terminated at this point. In Figure 39, the stress-strain plot of the 
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FIGURE 29. Contours of Constant Octahedral Shear Stress, Normalized by Dividing by the 
Matrix Yield Value of 17 ksi. Average Applied Stress , a * 97.4 ksi, 45* 
Section Longitudinal Model. x 
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flGUKu 30. Contours of Constant Octahedral Shear Strain (10 in/in). Average Applied 
Stress, a = 97.4 ksi, 45° Section Longitudinal Model. 




FIGURE 32. Contours of Constant In-Plane Shear Stress (ksi). Average Applied Stress, 
o =97.4 ksi, 45° Section Longitudinal Model. 
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FIMRE 36. Contours of Constant Maximum Principal Stress (ksi). Average Applied 
Stress , a - 183 ksi, 45° Section Longitudinal Model 
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DISTANCE IN FIBER DIAMETERS FROM PLANE OF FIBER BREAK 

FIGURE 30. Fiber Loadings versus Distance from Plane of Fiber BreaV, Average Applied Stress 
“ = 1 H 3 ksi 45° Section Longitudinal Model. 


AVERAGE APPLIED STRESS, a (ksi) 



FIGURE 39. Plot of Applied Stress versus Composite Strain for the 45° Section 
Longitudinal Model, 12.5 Percent Discontinuous Fibers. 



local composite Is presented, reflecting the occurrences of crack propaga- 
tion and Che associated decreases In composite modulus. 

The process of modeling stress redistributions due to crack growth 
involves many cycles of the numerical solution procedure, and when this 
1b coupled to the Incremental loading of an analytical model which experi- 
ences localized plastic deformation at a small fraction of its total load 
capucity, computer time requirements become very large. For example, the 
analysis of the 45° longitudinal section model discussed in this section 
required 9,577 seconds of running time on the CDC Cyber 760 computer and 
372 thousand words of core storage. When one considers that, in the 
present example, plastic deformation began at an applied stress level 
of 18 ksi, and at o ■ 208.5 ksi total failure had not yet occurred, 
the complexity of the problem being analyzed can be appreciated. 

6 . A . 2 . Crack Initiation and Propagation in the 90° Section Longitudinal 
Model 

The 90° section longitudinal model was loaded to a maximum of 208.3 
ksi, during which five periods of crack growth were observed, but catas- 
trophic failure of the composite did not occur. In the discussion that 
follows, contour plots and fiber loading diagrams are presented to illus- 
trate some of the more interesting points in the loading history of this 
model. While a great many similarities between the response of the 45° 
and 90° section longitudinal models exist, a careful comparison will in- 
dicate the differences in crack growth pattern, stress distributions, 
fiber loading, and load levels at which they occur. Such a comparison 
confirms the fact that there is a strong dependence of the results on 
fiber spacing when plastic deformation and crack growth are considered. 
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This dependence le what makes the use of a tvo-dinune tonal formulation so 
difficult to correlate to actual conditions, as there is no practical way 
of manufacturing most unidirectional composites in which every fiber le 
equidistant from each of its nearest neighbors. As was discussed in 
Section 4.1, it is possible that specialized test speclmena, used in 
conjunction with the axlsymmetric analysis, could lead to the determination 
of an average, i.e., effective, fiber spacing, thus allowing the present 
generalized plane strain formulation to more accurately model square and 
rectangular fiber arrays. 

Because of its much higher apparent fiber volume, shear loading of the 
matrix between the broken fiber and its nearest neighbor in the 90* 
section longitudinal model is more severe than in the 45° section longi- 
tudinal model. As a result, the initial stress concentration caused by 
the broken fiber is more pronounced, and initial element failure in the 
90° section model was observed before a plastic region of any significant 
extent could form around the crack tip. This failure occurred at an 
average applied stress level of 24.9 ksi. The matrix stress contour 
plots for the composite just prior to this event are presented in 
Figures 40 through 43. The single element failure at this applied stress 
level appears to have very effectively "blunted" the tu as the next 
failure did not occur until an average applied ? tress xevel of 84.0 ksi 
had been attained. At this point, one addition ' .lent failed, again 
temporarily blunting the crack. Figures 44 througn 48 illustrate the 
state of stress in the matrix and the loading of the various fibers at 
this load level just prior to the failure o f the * cond element. At this 
point there is a fairly extensive zone of jtic deformation around the 
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FIGURE 47 . Contours of Constant In-Plane Shear Stress (ksi) , Average Applied 
Stress ,o * 84 ksi, 90° Section Longitudinal Model. 
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crack tip. Also, the shear stress level* in the matrix between "liter 
fibers are much higher, and display some gradation. Of particular interest 
is the fiber loading nioc, Figure 48. In this 90* section model, the 
effect of the broken fiber is confined to a smaller region than for the 
45° section model (se Figure 28), with the adjacent fiber being loaded 
more than 20 percent higher at the plane of the break than the far field 
average. The more remote fibers are much less affected. At 8.2 fiber 
diameters, the broken fiber is almost fully loaded. 

The next stage of crack growth occurred at o ■ 85.4 kai, involving 
102 elements which failed in 9 intervals of constant stress crack 
propagation. The tip of the resulting crack advanced along the boundary 
of the adjacent fiber, as is also suggested by the final shape of the 
crack, shown in Figures 49 through 53. It will be noted in Figure 49 
thut the matrix has been unloaded to the point that once again there is 
no plastic deformation. In all these respects, the results of the 90° 
section model differ from those of the 45° section model. 

Loading was then increased monotonically to a level of 145.4 ksi, 
at which point one more element failed. Again, the states of stress 
and strain in the matrix and fibers have been plotted so that one can 
contrast them to the situation that existed just after this particular 
crack was formed. These results are presented In Figures 54 through 57. 

At thir, point there is a very extensive zone of plastic deformation, and 
the effects of the crack on the fiber loading are not nearly so localized. 

At a load level of 146.8 ksi, extensive crack growth again occurred. 

A total of 83 elements and 10 intervals of constant stress crack propaga- 
tion were involved in this process. The shape of the resulting crack, 
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FIGURE 49. Contours of Constant Octahedral Shear Stres Normalized by Dividing by the 
Matrix Yield Value of 17 ksi. Average Applied Stress, * 85.4 ksi, 90° 
Section Longitudinal Model. 
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FIGURE 51. Contours of Constant Maxinun Principal Stress (ksi), Average Applied 
Stress , 0 =85.4 ksi, 90° Section Longitudinal Model. 
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FIGURE 52. Contours of Constant In-Plane Shear Stress (ksi) , Average Applied Stress 
o = 85.4 ksi, 90° Section Longitudinal Model. 
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FIGURE 54. Contours of Constant Octahedral Shear Stress, Normalized by Dividing by the 
Matrix Yield Value of 17 ksi. Average Applied Stress,^ * 145.4 ksi, 90° 
Section Longitudinal Model. 
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FIGURE 56. Contours of Constant Maximum Principal Stress (ksi) , Average Applied Stress 
o = 145.4 ksi, 90° Section Longitudinal Model. 
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along with contours of normalized octahedral shear stress, are presented 
In Figure 58. It will be noted that once again the crack has grown through 
the previous region of plastic deformation, and little or no plastic 
deformation remains after this process. In Figure 59, the affect of the 
crack growth on the fiber loading can be seen. The two neurest fibers 
become affected near the plane of the fiber break, and in addition, the 
two more remote fibers begin to display the same pattern as the closer 
fibers, i.e., increasing load in the vicinity of the flaw, decreasing 
as the axial distance from the break site increases. 

Continued loading of the 90° section model to 208.3 ksi produced no 
additional element failures, at which point this computer run was 
terminated. The last states of stress to be plotted were at an average 
applied stress of 193.3 ksi. These are presented in Figures 60 through 
63. A plot of the composite stress versus strain for this example is 
presented in Figure 64. 

6.5. Crack Initiation and Propagation in the Transverse Section Model, 
Loaded T r ansversely 

Stress increments a were applied to the right hand boundary of the 
transverse model shown in Figure 11. Plastic deformation of the aluminum 
matrix began at the fiber-matrix interface at a point about 30° from the 
positive x-axis, i.e.. Element Nos. 88 through 91 (see Figure 11). As 
loading Increased, these elements became even more highly stressed. 

However, at 58.4 ksi. Element No. 81 ruptured in a hydrostatic tension 
failure node. This failure triggered the failure of Element Nos. 83 and 
85, whereupon further strain energy redistribution was required. Crack 
propagation progressed along the fiber-matrix interface, terminating at 
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FIGURE 61. Contours of Constant Octahedral Shear Strain (10 in/in). Average Applied 
Stress, a = 193.3 ksi, 90° Section Longitudinal Model. 
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the point shown In Figure 65. This failure pattern agrees very well with 
that presented by Adams [4] in his 1973 study of crack propagation in a 
transverse section of a unidirectional boron/alumlnum composite. His 
effort involved an analysis scheme with a constant displacement rather 
than u constant stress loading procedure, a finer mesh finite element 
representation, and smaller loading increments than were used in the present 
analysis. Adams also typically observed initial failure at matrix elements 
along the fiber-matrix interface at a point approximately 30° from the 
x-axls. However, his analysis did not Include a hydrostatic failure mode, 
which could explain the slight difference in the location of crack initia- 
tion observed in the present study. 

The loading increment following the 58.4 ksi load level resulted in 
further element failures, with the crack progressing along the fiber-matrix 
boundary. While this example was undertaken primarily for the purpose 
of comparing results of the present analysis with those of Adams [41, it is 
clear that the effects of disbonds, local matrix voids, and other manu- 
facturing flaws on the transverse strength of a unidirectional composite 
could also be characterized. With larger finite element arrays, constructed 
to simulate such flaws, manufacturing cycles of compression and thermal 
loading could be imposed to determine the residual stresses caused by 
these defects, and their effects on subsequent service loading and envi- 
ronmental exposure. 

6.6. Axial Loading of the Axisymmctric Longitudinal Model 

Numerous computer runs, involving a variety of finite element models, 
were made to verify the correctness and accuracy of the axisymmetric 


113 



FIGURE 65. Crack Propagation in the Transverse Section Model at an 
Average Applied Stress,? * 58.4 ksi. 
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analysis program under all possible combinations of element geometry, 

material types (i.e., elastic isotropic or transversely isotropic, and 

plastic, isotropic) and loading. When all difficulties had been identified 

and corrected, a series of runs was made, employing the finite element 

model shown in Figure 13. Three different variations of this model were 

employed in this preliminary study, aimed primarily at assessing the 

effect of thickness of the annular sheath of aluminum on the response of 

the single broken fiber configuration. The first model studied represents 

the condition of minimum fiber spacing in a 55 percent fiber volume, square 

array, unidirectional composite. This results in ;j model in which the 

broken fiber has a rather thin sheath of matrix around it, the ratio of the 

fiber radius to the total model radius (r c /r ) being 0.714. Another model 

r m 

studied was the case of maximum fiber spacing in the 55 percent fiber 
volume, square array composite, wherein the ratio of the fiber radius to 
the model radius is 0.424. Finally, after studying the results of runs 
using the first two models, xt was decided to Increase the thickness of 
the matrix in the maximum fiber spacing model by 50 percent, to further 
investigate the effect of matrix thickness in this particular configura- 
tion (r,/r * 0.330). 

r m 

The results of this study are summarized in contour plots for each 
model configuration, with normalized octahedral shear stress, octahedral 
shear strain, maximum principal stress, and in-plane shear stress being 
the parameters plotted. 

Figures 66 and 67 represent the case of minimum matrix thickness. 

It can be seen immediately that a rather small portion of the cross section 
of the model perpendicular to the applied load is aluminum matrix. Since 
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FIGURE 67. Contours of Constant Maximum Principal Stress and In-Plane 
Shear Stress, Minimum Fiber Spacing Model (r./r « 0.714), 

Average Applied Stress, o ■ 17.4 ksi. m 
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the average stress is applied to the entire circular cross section, there 
is an immediate stress concentration at the z ■ 0 plane (see Figure 13) 
due to an appreciable reduction in net cross-sectional area caused by the 
fiber discontinuity. This, and the presence of the penny-shaped crack, 
also a result of the fiber discontinuity, caused Initial plastic deforma- 
tion at the crack tip at an average applied stress of only 7 ksi. Continued 
loading resulted In a first failure in the matrix ai an applied stress 
of ]7.4 ksi. Failure of the first element triggered the failure of three 
additional elements along the z ra 0 boundary, the plane of the fiber dis- 
continuity. This resulted in a further reduction of net section at the 
z = 0 boundary, and a subsequent adjustment increment caused the crack to 
grow radially to the edge of the model, representing total failure and 
thus terminating the analysis. The contour plots of Figures 66 and 67 
represent the state of octahedral shear stress and strain in this model 
just prior to crack initiation. 

In the case of the axial loading of the maximum fiber spacing 
model (r^/r^ = 0.424), plastic deformation at the crack tip did not occur 
until an average applied stress of 13 ksi had been reached. Continued 
loading resulted in a first failure, or crack initiation, at a stress 
level of 32.3 ksi. As was the case with the minimum thickness model, the 
crack immediately grew through the region of plastic deformation, and 
the resulting reduction in the net section area of the model caused total 
failure. The contours of constant octahedral shear stress and strain 
for the maximum fiber spacing model just prior to crack initiation are 
presented in Figures 68 and 69. 
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OCTAHEDRAL 
SHEAR STRESS 



OCTAHEDRAL 
SHEAR STRAIN 



FIGURE 68. 


Contours of Constant Octahedral Shear Stress (Normalized by 
Dividing by the Matrix Yield Value of 17 ksi) and Octahedral 
Shear Strain, Maximum Fibei Spacing Model (r^/r m * 0.424), 
Average Applied Stress, o « 32 ksi. 
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For the thickest model run (r f /r u 0.330), plastic deformation was 

t m 

first observed at an averuge applied stress of 15 ksi. Initial crack 
formation occurred at an average applied stress of 37.2 ksi; the contour 
plots of stress and strain just prior to this crack Initiation are pre- 
sented in Figures 70 and 71. It win be noted that the region of plastic 
deformation around the crack tip involves more than a third of the matrix 
thickness, as was the case for the two models discussed above. As in the 
previous cases, this relatively low fiber volume model representation also 
suffered catastrophic failure once the crack had been initiated. 

To provide further insight to the response of the axisymmetrlc 
model, the displacement of the broken fiber ends relative to each other 
has been plotted in Figure 72 as a function of the average applied stress 
level for the three model configurations discussed above. This plot 
clearly illustrates the variation in overall stiffness and load carrying 
capacity of the model as a function of matrix thickness. Experimental 
measurements could readily be made using test specimens of these configu- 
rations, and results compared with those shown in Figure 72. 




FIGURE 71. Contours of Constant Maximum Principal Stress and In-Plane 
Shear Stress, 1*5 Maximum Fiber Spacing Model (r f /r « 0.33i 
Average Applied Stress , 0 "37.2 ksi. f ra 
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SECTION 7 

CONCLUSIONS AND FUTURE WORK 

AL the end of the second year of this continuing investigation for 
NASA-Lewis, the analytical tools required to perform the title study are 
now well in hand. Although Improvements can, and will, still be made in 
both the generalized plane strain and axlaymmetrlc analyses and related 
computer programs, they are presently fully operational. In addition, 
the three-dimensional finite element analysis, although developed as part 
of another program [7], has also just become operational. This analysis 
will also be fully available to the present NASA-Lewis study, as required. 
Obviously, analysis methods have advanced significantly during the past 
two years. 

The numerical examples presented in this report are intended to 
demonstrate the capabilities of the analyses, and to provide at least a 
preliminary indication of the influence of a broken fiber on local stress 
states, and overall performance of the composite. 

To date, only pre-existing fiber breaks have been modeled. It will 
be relatively straightforward to extend this to the analysis of composites 
having weak sites distributed arbitrarily along the fibers. The numbers 
and severity of these weak sites can be established from available experi- 
mental data for boron fibers. This will lead to the study of interactions 
between closely spaced fiber breaks occurring under an applied stress. 

In terms of maximizing energy absorption during plastic deformation 
and subsequent crack propagation, it is anticipated that a detailed sLudy 
using the generalized p’ 'e strain analysis will lead to guidelines for 
designing boron/aluir composites with controlled defects fabricated 
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Into the material. Thia will result In trade-off a between composite stiff- 
ness and ultimate strength, and energy absorption during the fracture 
process. The three-dimensional analysis can be used to confirm the ade- 
quacy of the (two-dimensional) generalised plane strain analysis In per- 
forming these studies. 

Experimental verification of the analytical work already completed, 
and to be undertaken during the next year, Is aico required. The axi- 
aymmetrlc analysis was developed primarily with this in mind. It will be 
relatively simple to fabricate single-fiber composites, l.e., single 
boron fibers surrounded by a uniform annular sheath of aluminum matrix 
material. Either a break can then be induced in the fiber before mechanical 
testing, or fibers known to have statistically weak sites can be used. 

When these single fiber composites are loaded in axial tension, the change 
in the gap between ends of the .*iber break can be experimentally monitored 
(using X-rays, an extensometer , etc.). As indicated in Figure 72, these 
changes are predicted to be relatively large, on the order of 0.002" to 
0.010" at failure. 

It was noted in Section 6.6 that crack initiation in the single fiber 
axisymmetric models lead to immediate catastrophic failure of the single- 
fiber composite. This was In contrast to the results presented in Section 
6.4 for the axially loaded, longitudinal section models using the gen- 
eralized plane strain analysis. There, the many surrounding unbroken 
fibers were able to absorb the energy released by the crack formation. 

Both analyses are presently set up to hold a constant average applied 
stress during the crack propagation and subsequent ‘djustment increments. 
Ba'ied upon the results of Section 6.6, it would be better to perform the 
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single-fiber composite experiments under displacement control rather than 
load control. Then, when plastic deformation, crack initiation, and 
subsequent crack propagation occurs, the average applied stress will drop, 
allowing the crack to be arrested. Crack opening displacement measurements 
can then be made before an additional increment of composite displacement 
is applied. 

Addition of a constant displacement loading scheme to the existing 
axisynimetric analysis will involve some modifications of the computer 
program. 

In conclusion, analysis methods are now well-established. Use of 
these analytical tools in performing detailed parametric studies remains 
to be completed. In conjunction with these analytical studies, experi- 
mental verification also remains as an important task. 
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APPENDIX A 


EVALUATION OF INTEGRAL COEFFICIENT S 
FOR THE GENE RA LIZED A XISYHMETRIC ELEMENT STIFFNESS MATRIX 

In Section 4.2.2, it was shown that a critical step in forming 
the exact element stiffness matrix for a toroidal finite element of tri- 
angular cross section is the integration of the product of the strain- 
displacement relationships and the constitutive relationships over the 
volume of the element, as described by Eq. (20). This operation led to 
the six integral relationships given in Eq. (24). In this appendix, a 
procedure for evaluating the three integrals having ^ terms in their 
integrands is presented. The three integrals to be evaluated are repeated 
here, i.c.. 


1 



drdz 


I 

I 


5 

6 



J J~ drdz 


(A-l) 


Tn Figure A-l, a planar section of 
reference to geometrical considerations 


the element is shown for 
in evaluating Equations (A-l). 
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FIGURE A-l. Geometric Definition of the Triangular Element in the 
r-z Plane. 

As Figure A-l indicates, the triangle, i, is defined by the lines 
bounding Its sides, i.e., L 12 , L 23 , and L 31 * each of which Is 
described by the equation shown. For example, line L 93 is expressed as 


2 “ ” 23 r + »„ 


(A-2) 


where m„„ is the slope of L _ and b-„ is its z-intercept, defined as 


. . Az . Z 3 " Z 2 

23 A r r 3 - 


< z 2 r 3 ~ z 3 r 2^ 
’23 (r, - r 2 ) 


(A-3) 
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The equations for the other two lines and their coefficients can 
be obtained by a cyclic permutation of the indices of Eqs. (A-2) and (A-3) , 
i.e., 1 -*• 2 •> 3 >1. Equations (A-l) can be directly Integrated with 
respect to z to yield 


I 


4 


I 


5 


I 


6 



(A— A) 


Equations (A-2) are substituted into Eqs. (A-4) to give 


f r 3i f V 2 r r 1 

r ( m 31 r+b 31> dr j r (m 23 r+b 23^ dr+ J < m 12 r+b 12 )dr 

f l r 3 r 2 


J 2r 

' i* 


<n. 3 i r +b 31 > dr+ 


(2 1 

| t lfe ( “23 +b 23 )2dr+ j r ^ (m 12 r+b 12 )Zd,; (A " 5> 


I “ 

e 


37< m 3i r+b 


3 l ) M r 2 i 


( r. 23 + b 2 3 ) dt+ 


A, 

J 3? (m 12 r+b 12 

^ r» 


) 3 dr 
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Equations (A-5) can now be expanded and Integrated term-by-term with 
respect to r fairly easily. The details of this operation are, however, 
quite lengthy, and the results only are presented below. In each case, 
only the first term is shown, the second two terms being obtained by a 
cyclic permutation of the Indices, as described In Section 4.2.2. 


h "I :b 12- b 31 )lnr l +<n 12'"31 ) 


I, - £ r< b 12- b 31 ) lnr l +(m 12 b 12- m 31 b 31 ) r l + 4 ( ”l2- B 31 )r l <**« 


3 <b X2‘ b 3X >lnr l + ' m 12 b 12' ,n 31 b 31 )r l + 2 <m X2 b X2" m 3X b 3X )r X 
+ 9 ( “X2" m 3X )r X + 


The expressions for 1^, I,., 1^ derived above are valid for the most 
general triangular geometries, i.e., when r^, rj, and r^ are distinct 
and not equal to zero. For certain orientations of the element, Eqs. (A— 6) 
are not valid; these orientations are described and dealt with in Appendix 
B. 

By substituting the expressions for the Blopes and z-intercepts of 
the various triangle sides into Eqs. (A-6) , then collecting and rearrang- 
ing terms somewhat, the expressions in Eqs. (26) through (28) are obtained. 
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However, It was found that computer programming of these Integrals Is 
easier and more direct when they are left In the form shown above. This 
la especially true when logical decisions regarding special element ge- 
ometry are programmed. 


APPENDIX B 


NUMERICAL DIFFICULTIES WITH SPECIAL GEOMETRICAL 
CONFIGURATIONS OF THE AXISYMMETRIC ELEMENT 


The expression for the six Integrals required for a full axlsymmetrlc 
element stiffness matrix, as given In Eqs. (25) through (28), are valid 
for the most general geometries, i.e., when r^, * and r^ are distinct 

and nonzero. However, there are three situations which require special 
treatment. These are: 

• When one of the node points of the triangle lies on the axis 
of rotation, i.e., r ■ 0. 

• When any two of the node point radii of the triangle are equal 
but nonzero. 

• When any two node point radii are zero. 

Each of these three cases are dealt with in the subsections that follow. 

B.l. One Note Point Radius Equal to Zero 

In examining the expressions for 1^, I,., and 1^ (Section 4.2.2 or 
Appendix A), it can be seen that logarithmic terms are involved. When 
a node point is located on the axis of rotation, its corresponding log- 
arithmic term becomes infinite. However, by examining the logarithmic 
term of interest and its coefficient, which consists of the z-intercept 
terms of the two lines converging on the node point in question, it is 
obvious that the intercept for both of these lines is the same. In other 
words, the limit of the logarithmic term and its coefficient can be shown 
to exist by an application of L'Hospital’s rule [17], and this limit is 
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always equal to zero. Accordingly* whenever the argument of a logarithmic 
term la equal to zero, the term and its coefficient are simply deleted. 

B . 2 . Two Node Point Radii Equal but Nonzero 

When two nodal radii of an element are equal, the element side 
between them la parallel to the axis of rotation, and the expressions for 
m^ and b^j, the slope and z-intercept of that side, become infinite. For 
example, if ^ ■ r^, an examination of Eqs. (A-3) quickly confirm the 
problem. This singularity is easily removed by considering the form of 
the Integrals in Eqs. (A-5), in which the integrations with respect to r 
have yet to be performed. Note that the integrals involving and b 2 j* 
the terms in question, also have limits of and r^. Thus, the integral 
is identically equal to zero and these terms can be omitted. In the 
implementation of the formulas given by Eqs. (A-6) , this objective is 
accomplished by defining the m and b terms of element sides parallel, 
to the z-nxis to be zero, i.e.. 


For r t - Ty 


m lj ■ b u * 0 


CB-l) 


B.3. Two Node Point Radii Equal to Zero 

When two of the element node points, say Node 1 and Node 3, lie on 
the axis of rotation, as indicated in Figure B-l, we observe 


r l “ r 3 " 0 

f -l - h ‘ 0 


U 1 ” U 3 


(B— 2) 


0 
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FIGURE B-l. Element Node Point Identification. 


When the above relations are substitute into Eq. (13) we have 
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(B-3) 


or 


t6 ^i " ^ T 00^i* 


(B-A) 
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and the assumed displacement field take* the form 
u - f, 2 r 

w - C 4 + S 5 r + £ 6 z (B-5) 

Applying the strain-displacement relationships of Bqs. (3) yields 

(B-6) 
or 

{l} i " tC 00 ]{f ' } i 

where [C Q0 ] is the shape matrix for an element with two node points whose 
radii are zero. By substituting this shape matrix into Eq. (21) and 
performing the indicated multiplication and integration, we obtain the 
element stiffness matrix for this geometry. In the case of an elastic, 
isotropic material, we have 
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Similarly, 


for a transversely isotropic material 


in the elastic range , 




. HI 
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0 I 0 I 0 I 0 

I JU-T+fMj I 0 I 0 

— r n 0“ 

1 o' 

intormc 


0 | 9 

o | <***r)t l 

~o 

0 I 0 

mfer'i i 0 

T (in)^ 


(B-8) 


whore the quantities Q, T, and F are as defined in Gq. (29). For an 
isotropic material in the plastic range we have, for two nodal radii 
equal to zero, 
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SYMMmiC 


4 ■ T>‘l 






(B-9) 


Where A, A', and B are defined in Eq. (30). 

The special form of the strain-displacement matrix, as given in 
Eq. (B-6), also requires that the back substitution matrices for each 
type of material response be re-derived. These are presented below. 
For an isotropic material in the elastic range we have. 
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00 1 


(1 + v)U - 


0 10 0 0 

0 2v 0 0 0 

0 0 0 0 (—p) 

0 v 0 0 0 


For the elastic, transversely isotropic case, 


I B oo’l 


“o (1-T+ F) 0 0 

E 0 <v+T+ F) 0 0 

Q 

0 0 0 0 

0 (1-T+ F) 0 0 


(v+T) 

(1-T) 


2(l+v) 


where T, F, and Q are as defined in Eq. (29). 

For an isotropic material in the plastic range, 
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with A, A', B, and having the same definitions as presented in 
Eq. (30). 



APPENDIX C 


LOAD APPLICATION IN THE DISPLACEK.^1 FORMULATION 
OF THE FINITE ELEMENT ANALYSIS 

A detailed description of the axlsyiranetrlc analysis was presented 
in Section 4. The generalized plane strain analysis was included in the 
first-year report [l]. Details of the computer programming were presented 
in Appendix A of that report. Of particular interest, however, is the 
method of load application. 

The present finite element analyses are displacement formulations, 
which is ideal in terms of accounting for the symmetry boundary conditions 
associated with the periodic arrays assumed. However, this presents a 
difficulty in terms of load application, since it is desired to be able 
to specify applied stress increments rather than applied displacement 
increments. In early works [4-6, 8, 9] this problem was handled by solving 
a series of displacement boundary value problems for each increment, one 
displacement boundary value problem for each component of loading incre- 
ment to be applied, viz, o , o , o , AT, AM. These individual solutions 

x y z 

were then scaled as required and superimposed to obtain the actual solution 
for the increment. Since it is the matrix inversion associated with the 
solution of each boundary value problem which consumes most of the computer 
time, doing this a number of times within each increment was very inefficient. 

Using a method introduced by Branca [10], it is possible to solve 

for any combination of mechanical (o , a , a ) and hygrothermal (AT, AM) 

x y z 

loadings in one step. This technique was incorporated into the basic 
raicromechanics analysis when it was first formulated by Miller and Adams [2], 




It continues to be used in the present program versions described in this 
report, having been refined and improved a number of times. A brief 
description of tbit "Branca" technique, as used in the present program, 
will be included here, for reference. 

The application of mechanical tractions to the finite element model 
Is considerably simplified by taking advantage of the rearrangement of the 
global stiffness matrix [K] , and the total force vector, { F } , using the 
method of Branca [10]. The displacement boundary conditions for the re- 
peating unit finite element model were specified in order to maintain 
continuity of the material continuum while satisfying symmetry requirements. 
Specifically, referring for example to Figure 6, displacements in the x- 
direction of node points along the right-hand vertical boundary must be 
uniform. Likewise, displacements in the y-direction of the upper hori- 
zontal boundary must be uniform, and the displacements of all node points 
in the z-axis direction must be uniform (the generalized plane strain 
condition). Whop the overall force-displacement equation of the system is 
considered, l.c., 

IF) - [K] {tf} ( C— 1 ) 

one can see that all of the boundary node points involved in mechanical 
loading will have identical displacements with respect to the direction 
of the load application. These identical displacements allow combining of 
certain terms in the global stiffness matrix that result in the replace- 
ment of the applied forces on boundary nodes by zeroes, in the manner 
described by Branca [10]. Successive modification of the global stiffness 
matrix for each boundary node point displacement results in the following 
form of Eq. (C-l) for the simultaneous application of uniform values of 
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where F , F , and F t are the total applied loads in the x, v, and z direc- 
x y 

tions, and are defined, for a unit thickness model, as 

F ■ a b 
x x 

F - o a (C-3) 

y y 

F = a ab 
z z 

where a and b are the lengths of the region of analysis (o.g., Figure 6) 
in the x and y directions, respectively. Modification of the global stiff- 
ness matrix is accomplished by summing the stiffness coefficients of un- 
known but equal boundary displacements throughout the system of equations. 
This results in a set of three equations, representing the forces along the 
three moving boundaries of the model, which are placed in the last three 
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columns of the global stiffness matrix, outside of the bandwidth of a 
normal stiffness matrix. Another summation is now performed on the coef- 
ficients in these three equations corresponding to equation numbers repre- 
senting loaded boundary nodes. These coefficients are all added to the 
terms in the last three rows of the three outside columns, which nov 
represent the total applied loads in the x, y, and z directions. The 
system of linear, simultaneous equations that results from this process 
involves a stiffness matrix that is no longer symmetric, and whose 
bandwidth has been violated. To further complicate matters, the banded 
portion of this stiffness matrix must be stored in rectangular form, as 
described by Zicnkiewlcz [12], to minimize the core storage requirements 
of the system of equations. The coefficients for the summit ten of force- 
equations are stored along side the rectangularized , upper triangular 
portion of the stiffness matrix. 

This system of equations is solved using a highly specialized form 
of Gaussian elimination in which the stiffness matrix and the load vector 
must be further modified. This solution technique requirec a great deal 
more bookkeeping than is the case for more conventional applications of 
the Gaussian elimination technique. The important advantage of this 
procedure, however, is that it allows the simultaneous application of 
external tractions in all three coordinate directions, and the application 
of thermal or moisture loads in a single step. 



